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In this supplementary paper we present some details on the solid-liquid interface detection, the 
deduction of the non-equilibrium free energy, the analysis of the granular temperature and energy 
per mode, a validation of the small slope approximation, a description of the Langevin dynamics, 
and the error analysis. 

A. Experimental determination of solid-liquid interface 



o 

; i In the quasi-2D geometry the solid phase consists of two square interlaced layers instead of the single hexagonal 

Qh layer that is characteristic of 2D systems [1]. These interlaced square lattices, when projected in 2D, result also in a 

square lattice. If particles in the solid phase are permanently in contact such that grains are close packed, then the 

unit cell per layer should have an area d x d, implying that the projected lattice should have a Voronoi area equal to 

d 2 /2. 

Two local parameters can be used to analyze the system. For each particle we can compute its Voronoi cell area, 

A v , which is inversely proportional to the local density, and also the 4-fold local order parameter, Q 4 , that is defined 
O as [2] 

*^J Here, Nj is the number of nearest neighbors of particle j and a{ is the angle between the neighbor s of particle j and 
the x axis. For a particle in a square lattice, |<3 4 | — 1 and the complex phase measures the square lattice orientation 
respect to the x axis. 

Fig. IT] presents the probability distribution function (PDF) of normalized Voronoi areas 2A V /d 2 and local order 
parameter |Q 4 | obtained from 3000 images, for which the system is indeed phase separated. The Voronoi area, A v , is 
normalized by d 2 /2, the corresponding Voronoi area for closed packed two square interlaced layers. Fig. ^1 shows the 
marginal PDF, defined as 



PDF(A V ) = JpDF(\Q 4 \,A v )dA v , (2) 
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(N 

Tfr PDF(|Q 4 |) = /PDF(|Q 4 |,A„)d|Q 4 |. (3) 

o J 

From the joint and marginal PDFs (Fig. [Tland Fig. [2] respectively) we can state that there is clearly a superposition 
of three distributions: a wide distribution related to the liquid phase, around |Q 4 | « 0.2 and 2A v /d 2 « 2; another 
more localized and stronger peak related to the solid phase, close to IQ4I = 1 and 2A v /d 2 « 1.3; and a third small 
peak of more ordered particles in a dense phase, with |Q 4 | w 0.55 and 2A v /d 2 sa 1.4, which correspond to small 
density and order fluctuations in the system. In general, the normalized Voronoi area of particles in the solid phase 
is larger than 1, implying that particles are slightly separated (in average). A separation distance of 0.15d is enough 
to shift this peak from 1 to 1.3. 

For classifying particles in the solid or liquid phase three simple criteria can be established. The first one is to 
use only A v , defining a critical value A c v (2A c v jd 2 « 1.5) such that particles with A v lower (larger) than this value 
are labelled as solid (liquid) particles. A second criterion is to use only local order. If |Q 4 | ^ Q\ (|Q 4 | < Q 4 ), then 
the particle is in the solid (liquid) phase. A third possibility is to use both quantities, such that if A v ^ A% and 
IQ4I ^ Q% are both satisfied, then the particle is considered to be in the solid phase. We have verified that our results 
are not sensitive to which condition we use (see section [M), and for simplicity we choose the second condition, with 
\Q%\ = 0.7. However, the first criterion also detects particles that are less ordered but dense, adding a thin liquid-like 
layer, with intermediate order, around the solid cluster. 

We note also that in Ref. [2] we concluded that Q 4 is a more appropriate order parameter to characterize the liquid 
to solid transition, compared to the local density. Indeed the correlations of Q 4 showed critical divergencies at the 
critical point, while those of the local density were insensitive to the transition. 
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FIG. 1: Probability density function of normalized Voronoi areas 2A v /d 2 and local order parameter IQ4I obtained from 3000 
images (about 2.8 x 10 6 pairs of A v and IQ4I values), for which the system is indeed phase separated. The criterion for solid- 
liquid identification relies on the valley obtained around JQ4J = 0.7 that is almost independent of A v . The color code indicated 
in the vertical bar is the logarithm (base 10) of the PDF. Contour plots are shown for log 10 [PDF(|Q4|, A v )] = —2, — 1, —0.5, 
0, 0.1, 0.2, 0.3, 0.5, 1. 
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FIG. 2: Projected probability density functions of normalized Voronoi areas 2A v /d? (a) and local order parameter IQ4I (b) 
obtained from 3000 images. The minimum obtained for IQ4I = 0.7 is chosen as criterion for solid-liquid identification. 



Figure 1 of the main text shows a typical example of the solid-liquid interface detection. In our experiment, for 
each image particles that are in a solid or liquid phase are labeled using the previously defined criterion. Also for 
each image, particles belonging to the same cluster are identified, as well as the largest cluster in each image. This 
cluster of course corresponds to the largest, most stable, solid cluster that is present when the driving acceleration 
exceeds the critical value. Then, for each image the center of mass of particles that belong to the largest solid cluster 
is determined. Using all the centers of masses from a stack of images, the global center of mass is determined. 



The next step then corresponds to detect those particles that are at the solid-liquid interface of the largest solid 
cluster for each image. This is done by an angle coarse-graining procedure. From the position of the global center 
of mass, the position of the solid particle that is the most distant from this center of mass within an angle A9 is 
determined. This is done from $i = 0° until 6 N = 360° - A6, following 9 t = 0i + (i - 1)A0. Thus, for each i; there 
is a pair of coordinates {x\,y\) of the particle that is in the largest solid cluster and that is the most distant to the 
global center of mass, thus at the solid-liquid boundary. For simplicity, we describe each of these particles by the 
coordinates {0\,r\), which are polar coordinates of each particle respect to the reference system fixed at the global 
center of mass. The final step is to interpolate these coordinates to a set where the angles are equally spaced, such 
that 9\ nt = 6\ + (i — l)A0/2, and we end with the set of coordinates (0™*,)-™*). It is this set of coordinates that is 
used for the Fourier analysis presented in this study. 

Figure 1 of the main text also shows that in general the interface is rather smooth, except for some "jumps" . 
These jumps are artificial, as they arise because we force the interface to be described by a single-evaluated function 
r l nt — r i nt (^ nt ); which is not always the case. Because of the adopted procedure for the interface determination, 
we necessarily end with a few of these jumps for each image. However, theses jumps do not contribute to the low 
wavenumber modes that turn out to be relevant for this study. This is explicitly shown in section [FJ where we discuss 
the validity of the small slope approximation. 

B. Non-equilibrium free energy 

The solid-liquid interface non-equilibrium free energy of a curved surface is written in two dimensions as 

E = 7 J ^B 2 + (d e R) 2 d9. (4) 



The total non-equilibrium free energy has an additional term, because if we minimize the solid-liquid interface non- 
equilibrium free energy, the absolute minimum is at R = 0, that is, no droplet. In principle, this is solved by adding a 
Lagrange multiplier that fixes (in average) the total mass (or equivalently, the area). Then, the total non-equilibrium 
free energy should be 

2w 2tt 

E = 1 [^/R 2 + (d e R) 2 d6 -nj^-dd. (5) 



where fj, is the Lagrange multiplier, which apart from a proportional factor is equivalent to the effective chemical 
potential. The negative sign is arbitrary and it is there to further simplify the notation. 

The equilibrium droplet is obtained minimizing the non-equilibrium free energy. First, non- uniformities in the 
radius increase the non-equilibrium free energy, so we assume a circular shape R = Rq. Then the non-equilibrium 
free energy is 

E = 27T7i? - HirRl- (6) 

To extreme E we take the derivative and equal it to zero. The solution is 

Ro = 1 , (7) 

or equivalently /x = j/Rq. Note that this last expression shows that we have correctly chosen the sign for /i in fcfy. 

Now, we linearize about Rq. That is, R{&) = Rq + eSR(6) and keep up to quadratic terms in e in the energy. Doing 
the expansion we obtain 

2-k 2-k 

E = i (\r + €5R + e 2 iM^ 1 d9 - g I [Rl + 2eR SR + e 2 SR 2 ] dS. (8) 

J 2/t 2 J 



Recalling that /1 = j/Rq the linear terms cancel (as it always should happen when doing an expansion about an 
equilibrium state) and the final result is 

2ir 

E = 7r 7 i?o + ^r f [(d e 6R) 2 - SR 2 ] d9 (9) 

^Rq J 




This expression does not describe our system. Indeed, the non-equilibrium free energy increases due to deformations 
of the interface but is decreases when changing globally the radius. That is, i?o is unstable under changes of the 
radius. In fact, this inconsistency could have been predicted when looking at the expression Q: it is clear that R is 
a maximum and, therefore, un unstable equilibrium. 

1. The origin of the problem 

The problem originates in the election of the thermodynamic ensemble. Indeed, when using fi we moved to the 
grand canonical ensemble in which the number of particles is not fixed. What we have found is what is called the 
Critical Nucleus in the Homogeneous Nucleation Theory. This is the size we have to overcome to create a droplet 
that, after reaching that size, will grow indefinitely. Droplets smaller than this radius will shrink again due to the 
energy cost of the free surface. The growth without bound is possible because we are in the grand canonical ensemble 
where we have fixed // and, therefore, there are always available particles to change phase. 

2. The solution 

The solution consists on changing the ensemble to one in which we fix the total number of particles. This is a 
complicated ensemble to make calculations because we have to take into account the physical fact that when particles 
are moved to one phase to another, the supers aturation changes and therefore the chemical potential is dynamically 
adjusted. This is precisely what happens when, working in the canonical ensemble, we prepare the system with a 
density between the liquid and solid densities. Clusters of the solid phase will be created, decreasing the density of 
the remainder part until it reaches the liquid density. 

To model this solution, we can write in the case of a circular droplet 

E = 27r7.R0 - fi{Ro)irRl. (10) 

If /j, is a decreasing function of i?o, a second equilibrium appears (this time stable) at a larger value of Rq- This is the 
final equilibrium radius of the droplet. 

C. Full model 

Equivalently of using a chemical potential that depends on the radius, we modify (I5J) to have a non-linear dependence 
on the total mass of the cluster, expression that will be easier to manipulate 



2tt / 2-7T \ 

E = 7 y Vi?2 + {d 6 R)H6 -flj ^dO j 
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(11) 



where / is a nonlinear function. For a cluster of area much smaller than the total system area, f(x) — /ix as in (fsl) . 
Again, an expansion up to quadratic terms is done obtaining 



E = 7 

o 



ZHq 
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- (f {irRl) + ef (irRl) R h + - [/' (irl%) h + /" {nR 2 ) R 2 a h] J , (12) 
where 

h = f SRdO (13) 

h = ({SRfaie (14) 



/ ; = j / 8Rd0\ . (15) 



The equilibrium radius is such that the linear terms should cancel (as to have a global minimum). This condition 
gives 

f(irR 2 )R Q = 7, (16) 



that once substituted in ( 12 ) gives 



2tt / 2tt \ ^ 

E = Eo + ^J [(d e 8R) 2 ~ SR 2 ] d6 - ~f"(irR 2 )R 2 J 5Rd0 . (17) 



o 

Fourier-transforming the last expression gives 



where 



E = E + ^P- V \SR m \ 2 (m 2 - 1) - ^-r( 7 ri?2)(2^i? ) 2 | ( 5i?o| 2 , (18) 



SR m = ^ I SRe-" n9 d6 (].!)) 

Z7T 



In the expression (18) all terms can be correctly interpreted. First, each mode contributes to the non-equilibrium 
free energy with a curvature dependence to 2 , proportional to the surface tension. However, the modes m = ±1 do not 
contribute. This is indeed correct because when a circle is slightly (linear in e) perturbed in the modes to = ±1, it is 
only translated but not deformed (this is related to the absence of dipolar perturbation of a circle or sphere) and can 
be directly checked as follows. Take the radius and do a to = 1 deformation, for example R = Rq(1 + ecos(8)), then it 
is simple to check that (x — eRo) 2 + y 2 — R 2 , that is the circle is just translated. As there is translational symmetry 
these modes should not contribute to the energy. Finally, the last term corresponds to the increase of energy when the 
mass of the cluster is changed. To assure that the equilibrium is stable, the second derivative of / should be negative 
and its absolute value larger than the first term. Note, however, that its value is otherwise completely independent 
of 7, that is the mass fluctuations are independent of the shape fluctuations. In summary, the non-equilibrium free 
energy about the equilibrium droplet can be written as 

E = E + ^\SR \ 2 + P V \6R m \ 2 {m 2 - 1), (20) 

R ° R ° |H>2 

where the formal expansion parameter e has been suppressed and A is a new parameter, which has the same units as 
7. In experiments, the translational symmetry is not perfect and the cluster has a tendency to remain in the center 
of the box. Then, the non-equilibrium free energy should be modified to 



E = E Q + -^\6R \ 2 + —{\6R 1 \ 2 + \SR, 1 \ 2 )+-L V \6R m \ 2 {m 2 -l) (21) 




with v a new parameter, with the same units as 7 and A. This expression is used to derive the equilibrium power 
spectrum and the time correlation functions. 

The parameters A and v have physical origins that arc different from the surface tension parameter 7, so their three 
values do not have to be directly related. A corresponds to a measure of solid cluster's size changes, either by expansion 
and contraction (at a fixed number of particles) or by the condensation and evaporation of particles in and out of the 
solid cluster, v corresponds to a measure of how far from the ideal condition is the current experimental realization. 
Ideally, v — 0, which implies that this mode (translation) is a neutral mode, equivalent to the spatial average for a 



flat interface (limit k — > 0). If v = 0, the m = 1 mode would realize a simple random walk, as discussed in section G 
In fact, we show below that this is not the case. We speculate that the origin of v > is related to particle's vertical 
dynamics, and thus depends on surface roughness, friction, and top and bottom local wall parallelism. In particular, 
it depends on the vertical dynamics of the whole solid cluster, which is still largely unknown. 



D. Static power spectrum 



Starting from (21), for equilibrium systems the equipartition theorem from statistical mechanics is invoked. Thus, 
each normal mode, when configurationally averaged (in our case when time averaged), contributes ksT/2 of energy 
to the system. Thus, the power spectrum is 

(l^ol 2 ) = k ~f^, (22) 

QSR ±1 f) = *^S, (23) 

<i«-' 2 > - S^iJ- < 24 > 

where the last line is valid for \m\ ^ 2. 

In general, the equipartition theorem is used for the total kinetic and potential energies (if any) by defining the 
number of active modes n. Thus, in average 

(K) = (U) = n h B T, (25) 

where K and U are the total kinetic and potential energies respectively. The idea is that active modes contribute 
with kgT/2 to the total energy, whereas non-active modes do not contribute. The number of active modes n depends 
on the particular system, by the number of spatial dimensions and by the particular atom or molecule interactions 
(intramolecular and intermolecular) . The most simple example is an ideal monoatomic gas in D s dimensions, with 
no potential energy. In this case, n = D S N , where N is the number of atoms. 

In our case, we consider that the granular system has an effective "granular" temperature T e ff. Although this 
quantity is commonly referred as a temperature it has units of energy. A simple version is to consider it equivalent 
to ksT. But, our non-equilibrium system shows that this effective thermal energy, or granular temperature, has to 
be defined carefully, as well as the number of active modes. In our setup, the system is isotropic in the horizontal 
plane, but the horizontal thermal energy is different from the vertical one. For the following analysis we focus our 
attention to the projected 2D dynamics, which is what can be analyzed experimentally. In fact, figure 2b of the paper 
shows that equipartition can be applied for the horizontal kinetic energy, up to a given wavenumber (from m — 2 to 
mw 30). 

We now focus on the solid-liquid interface, as a subsystem of the complete granular system. We recall that the 
protocol used identifies one interface particle for each angle interval A6 — 2°. Therefore, the interface subsystem 



consists always in exactly N p = 180 particles. The interface total potential energy is given by equation (21 ). We then 
consider equipartition for the total kinetic and potential energies, 

(K) = {E)=n l -T cS , (26) 

where n is the number of active modes and T B g is the effective granular temperature. The important issue is that we 
can measure (K) and express it as a sum of normal modes. Indeed, 



W = \ E m P^) = l N P m p(v 2 ) (27) 



i.=i 



where (v 2 ) — (v 2 ) + (v 2 ) is the velocity variance. Here, we must stress that because v is function of both 9 and t, 
the average ( ) is a double average, over angles and over time. Both velocity components can be expressed by Fourier 
series: 

oo oo 

«x= E r m e m<> , v y = J2 < eim<> > ( 28 ) 

m— — oo m— — oo 

where v^ and v^ are the Fourier components. In practice, because A9 = 2°, then the summation is done for a finite 
number of modes, from to = —to* to to = to*, with to* = N p /2 — 1 = 89. Thus, 

(K)= 1 2 N p m p J2 (Kl 2 + W)> ( 29 ) 



which allows to define the kinetic granular energy per mode, 

(Km) = \N p m p (\v x m \ 2 + Kf). (30) 

For the last two expressions, the angle average has already been performed, and ( ) is a configurational (time) average. 
This is the quantity plotted in figure 2b of the paper, which shows equipartition from m — 2 to m f=a 30 (corresponding 
to wavelengths rj 74d to w 5d). In fact, an average between m = 2 and m = 30 gives (K m ) = K oq — 4.82 ± 0.04 nJ. 
The two lowest modes, which correspond to global radius fluctuations and center of mass translations, have higher 
energy: (Kq) = 8.0^4'g nJ and (Ai) = 5.9 ± 2.8 nJ. Their errors are slightly asymmetric due to their asymmetric 
(non-Gaussian) distribution functions, which are detailed in the last section of this document. 
Furthermore, the identity 

1 1 m * 

(A) = -N p m p (v 2 ) = -N p m p £ (|<| 2 + ^| 2 ), (31) 



m— — m v 



can be verified numerically. Indeed, from the particle velocity distributions (in x and y) and their variances, we obtain 
(A') sa 0.80 /xJ, whereas from the velocity Fourier components we obtain (A) sa 0.82 /iJ. 

Next, as we measure for each mode (K m ), we can apply (K m ) = (SE m ), where (SE m ) corresponds to each mode 
contribution to 



(SE) = ^(\SR \ 2 ) + ^ ({\SR!\ 2 ) + (|^-i| 2 >) + 7 E d^™| 2 )(™ 2 " !)■ (32) 



|m|>2 

Thus, the final expressions for the static power spectrum for our system are 



(|«?o| 2 ) = (K °^°, (33) 

(\6R ±1 f) = ^1^, (34) 
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irj(m z — 1) 



i«w 2 i 2 ) = n rf ; ■ (35) 



E. Capillary-like spectrum for different interface detection conditions 

In section A we have presented the interface detection procedure, including the three different criteria used for the 
selection of particles in the solid cluster. In the main text we also discuss the validity of the capillary-like spectrum 
respect to the choice of the coarse-graining angle A8. 

In figure k^ we present several fluctuation spectra (|<5i? m | 2 )/(A m ) versus m 2 — 1 using the three different criteria 
and several coarse-graining angles. The l/(m 2 — 1) dependence is observed for all cases, with variations on the range 
of mode numbers m for which is it valid, depending on the specific details of the detection procedure. For Ad = 2° 
(Fig. [3k) we obtain that for the first criterion the fitted surface tension results 7 = 2.0 ± 0.1 /iN. For the second and 
third criteria, within errors their result are equal, 7 = 2.9 ± 0.1 /iN. Figs. [3]d and[3J: present spectra for the second 
criterion with different coarse-graining angles. 

These results allow us to conclude that the description of the fluctuations by the capillary theory is very robust 
with respect to the details of the interface detection procedure. However, the numerical parameters of the theory 
(surface tension and mobility) are more sensible to the different criteria. In particular, the fact that 7 is lower for 
criterion 1 can be explained recalling that the detection procedure adds a thin liquid-like skin to the solid cluster, as 



mentioned in section A This layer is composed by particles with intermediate order (IQ4I « 0.55) but with small 



Voronoi area (2j4„/cP ~ 1.4), thus with higher density. 

F. About the small slope approximation 

The small slope approximation is used explicitly in the Taylor expansion of \j ' R 2 + (dgR) 2 . In particular, the last 
two terms in 
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FIG. 3: (a) Fluctuation spectra {\8R m \ 2 )/{K m } versus m 2 — 1 for AG = 2° and the three cluster selection criteria presented 
in section fSl The vertical dashed lines show the limit values used for the fits: m — 6. ..30 for criterion 1, and m = 6. ..54 for 
criteria 2 and 3. (b) Fluctuation spectra {\SR rn \ 2 )/{K rn } versus m 2 — 1 for several A8 using criterion 2. For smaller A9, noise 
contaminates the spectrum at higher m values. For larger A8, the coarse-graining procedure acts as a filter at larger m. (c) 
Same as (b) for a smaller range of m, putting in evidence the discrepancies that appear for larger m. 
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(36) 

(37) 



are neglected, as well as the higher order terms (h.o.t) that are not shown. As discussed previously, the contribution 
of the linear term vanishes in the non-equilibrium free energy. Thus, we are left to show that the following conditions 
are indeed fulfilled by our experimentally determined solid- liquid interfaces: 
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(38) 
(39) 



These terms contribute to the non-equilibrium free energy through the angle integral. Thus, by expanding 5R in 
Fourier series, the above expressions become equivalent to show these conditions are satisfied: 
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(40) 

(41) 



Here, the brackets ( ) represent an ensemble average, which in practice is done through time average of data from 
many images. These conditions are indeed satisfied, as shown in Fig. [IJ The first parameter satisfies ^i < 1.7 x 10~ 6 
for all to, whereas the second one, ^ < 8 x 10~ 3 for all to. 
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FIG. 4: Relative intensity of the correction terms, *i = (|<5#T l | 4 )/(8-Ro) (left) and * 2 = m 2 (\5R^\ 2 ) / (ARq) (right), to the 
quadratic form of the non-equilibrium free energy in (37 1. Small values imply that the small slope approximation to the 



non-equilibrium free energy is valid and the quadratic form can be used. 



G. Dynamic correlation function: determination of effective mobility 



We now focus on the mobility parameter M and the interface dynamic correlation function. For an interface of flat 
geometry, with local height h(x,t), the standard procedure is to consider the interface local velocity V = dh/dt = 
Mjk, where k is the interface's curvature (/c — d xx h for a flat interface in 2D). The Fourier transform of this equation 
leads to 



dt 



= -M^k 2 A k , 



where we have used the Fourier representation 



h( X ,t) = j2Mty k 



(42) 
(43) 



When a noise term is included, Langevin equations are obtained, in real and Fourier space respectively: 



dh(x,t) / d 2 h(x,t) 



dt 



dx 2 



where 



dAk 
dt 



V (x,t) = J2m(t)e ikx 



M(--yk 2 A k + r) k (t)), 



(44) 
(45) 



is modeled as a white noise term, which is delta correlated 



( V (x,t)) = 0, 

(r](x,t)r](x' ,t')) = C5(x-x')S(t-t'), 

(Vk(t)) = 0, 

(Vk(t)r,* k {t')) = C6 k ,«6(t-t'). 



(46) 
(47) 
(48) 
(49) 



In 2D, if we consider an interface of length L between two semi-infinite phases, then the non-equilibrium free energy 

is 

(50) 

k 



E = j[ Vl + {d x K)Ux -lL+^Y] k 2 \A k \ 
Jo z i. 



It is straight forward to show that Eqn. ( 42 ) can be obtained by the general form 

dA k M SE 



dt 



L SA k 



10 



(51) 



where SE/5Ak is the functional derivative of E. 

Applying the same formalism to our case, and including noise, from Eqn. (21 1 we obtain the following Langevin 
equations 



dSR 
dt 

dSR ±1 



dt 



M 



M 



rj^5R + r) Q (t) ) , 



i' 



SR±i + v±i(t)) , 



dSR 



= M[-^(m z -l)SR m + r lm (t) 



"o 

dt v JL0 

These Langevin equations have the form x = —x/t + Mrj(t), which solution is 

x(t) = x{fd)e~ t ' T + M f e- {t ~ s)/T ri{s)ds. 
o 

The ensemble average results (x(t)) = (a;(0))e _t / T , and (x(t)) — > for t — > oo. 
The two time correlation function is 



(x(t)x*(t')) = (\x(0)\ 2 )e-( t+t 'V r + ^^(< • ■ ■ -<■ (M; ' 



I - e- 2t /-' 



For t = t' ', the static power spectrum is used in long time limit of 

2,-2t/r M 2 CT 



(\x(t)\ 2 ) = (\x(0)\ 2 )e-^ + 



which gives the constant C = 2(K m ) / \itRqM) . 
Finally, the mean square displacement is 



From Eqn. (55) we obtain 



(|A J R| 2 ) = (|x(t)-x(0)| 2 ). 



x{t) - x(0) = x(0)(e- t/T - 1) + M / e-^- s ^ T T]{s)ds, 



(52) 
(53) 
(54) 

(55) 



(56) 



(57) 



(58) 



(59) 



thus 



(|x(t)-x(0)| 2 ) = (|a;(0)| 2 )(e- t / r -l) 2 + M 2 / / C -(M— «')A ^(s)^ '(s')}dsds' 







(I^^Ce-^-l^ + ^^fl < 2 



(60) 



(61) 



where we have used the fact that {x{Q)rj(s)) — (x(0))(r](s)) = 0, for s > 0, and the double integral is solved as before. 
We remind that the long time limit of (57) gives M 2 Cr/2 = (|a;(t)| 2 ) = (|x(0)| 2 ), allowing to simplify expression (61 1, 



leading to the following expression for the mean square displacement 

(\x(t) - x{0)\ 2 ) = M 2 Ct (l - e- f/r ) . 
Finally, replacing the expression for C in terms of the energy per mode, 



(62) 



(63) 
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FIG. 5: (a) Probability density function, PDF, and (b) cumulative density function, CDF, of Kq for one realization. In (a), the 
continuous line shows a Boltzmann distribution fit, PDF(i\ ) = Aexp(—K /B), with A = 0.14±0.02 (nJ) -1 and B = 8.1±0.4 
nJ. In (b), the vertical solid line corresponds to the average (Ko) — 8.5 nJ; the vertical dashed lines show the limits for which 
16% of the data is below (Ko) — &-, and 16% is above (Ko) +cr+, with <r_ = 6.9 nJ and <r+ = 7.1 nJ. For the other realization, 
(Ko) = 7.5 nJ, a- — 6.2 nJ and <r+ = 6.1 nJ. 



H. Error analysis 



Both |<5i? m | 2 and K m have large fluctuations. Their distributions are strongly non-Gaussian, which implies that the 
corresponding error analysis has to be realized carefully. For example, the measured average kinetic energy for m = 
is 8.0 nJ, whereas its standard deviation is 7.9 nJ. Thus, the error related to this measurement is not its standard 
deviation. In fact, the associated error has to be computed with a generalized criterion. When a measurement of a 
quantity x is performed, the usual procedure is to assign the standard deviation <r x as error to the average (x) of a 
set of data. For a Gaussian distribution, the range [(x) — a x , (x) + a x ] corresponds to a confidence interval of ss 68%, 
meaning that if a new measurement is performed, it has a 68% probability to be in this range. 

Figure [5^, shows that Kq obeys a Boltzmann distribution, strongly non-Gaussian. Fig. [5]d presents the cumulative 
distribution function (CDF) of the data presented in Fig. [5k. The adopted criterium is the following: errors er_ and 
ct + arc computed by estimating the values of the measured quantity that insure that 16% of the data is lower that 
(Kq) — <t_, and 16% is larger that (Kq) + a+. Thus, 68% of the data lies in the range [(Kq) — <r_, (Kq) + a + }. The 
same procedure is realized for each mode m. 

Figure [6^, shows that K\ obeys a generalized Poisson distribution, also strongly non-Gaussian. We have not yet a 
satisfactory explanation of why this fit seems to work very well, so it must be considered as a phenomenological fit. 
The distribution's asymmetry affects less the measured error bars. For the particular realization shown in this figure, 
<t_ = 4.1 nJ and er + = 4.1 nJ, with (K\) = 6.1 nJ. However, the average is clearly larger than the most probable 
value K" lp = 3.0 nJ (also known as the maximum likelihood value). For m ^ 2, all K m obey the generalized Poisson 
distribution. Furthermore, between m = 2 and m s=s 30 they all collapse on a single curve, which is consistent with 
the observed equipartition. 

The fluctuations of |5i? m | 2 also obey similar non-Gaussian, asymmetric distributions. Their error bars are computed 
using the same procedure described above. 

For both |<5i? m | 2 and K m the final asymmetric error bars are computed by averaging two independent realizations, 
thus reducing the total error by a factor s=s v2. This is what is known as the standard error for N independent 
measurements, which scales as a/\/N, where a represents the error of each independent measurement. Because the 
asymmetries in the errors for both |W? m | 2 and K m are small, the reported error bars in the paper correspond to the 
maximum of <r_ and er + , which allow an easier computation of errors when quantities with uncertainty are combined 
(sum, difference, multiplication or division) or when they are used for computing another quantity by means of a 
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FIG. 6: (a) Probability density function, PDF, and (b) cumulative density function, CDF, of K\ for one realization. In 
(a), the continuous line shows a generalized Poisson distribution fit, log[PDF(A'i)/d] = ~a[K\ + b) + c\og[a(Ki + 6)], with 
a = 0.32±0.04 (nJ) -1 , b = -0.02±0.60 nJ, c = 1.0±0.5 and d = 0.33±0.06 (nJ)" 1 . In (b), the vertical solid line corresponds 
to (K\) = 6.1 nJ; the vertical dashed lines show (Ki) — <r_ and (Ki) + <j+, with <r_ = 4.1 nJ and a+ — 4.1 nJ. For the other 
realization, (Ki) — 5.7 nJ, <r_ = 3.8 nJ and a+ = 3.8 nJ. 



non-linear function. 
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